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Abstract 


This report discusses some analytical procedures to enhance the real time solutions of PMARC 
matrices applicable to the Wall Interference Correction Scheme (WICS) currently being implemented 
at the 12 toot Pressure TunneL WICS calculations involve solving large linear systems in a reasonably 
speedy manner necessitating exploring further improvement in solution time. This paper therefore 
presents some of the associated theory of the solution of linear systems. Then it discusses a 
geometrical interpretation of the residual correction schemes. Finally some results of the current 
investigation are presented. 


Introduction 


WICS is a combined experimental and computational approach to correct for the wall interference 
effects in wind tunnel testing. The procedure involves replacing the test model by a combination of 
mathematical singularities. The use of the model force and moment measurements together with the 
wind tunnel wall pressure measurements facilitate calculation of the unknown strengths of sources 
and doublets representing the model and the support system The main benefit of the WICS procedure 
is in the pare-calculation of a database using the PMARC influence coefficient matrices. Often this 
process takes several days depending upon the number of singularities and panels in the model. 
Moreover the influence coefficient matrices arising from WICS calculations are neither sparse nor 
symmetric. Thus exploring a better linear solver would be a worthwhile effort. With this objective in 
mind, this investigation attempted to answer some of the related issues of the solution of linear 
systems. 

This report re-visits some popular direct and iterative methods of solution of linear systems and 
systematically compares their performance for solution of large systems. For some time the 
convergence difficulties of PMARC matrices applicable to WICS calculation was believed to be ill- 
conditioning of the influence coefficient matrices. Thus the early part of this research was devoted 
to the study of ill-conditioning. As explained later, this issue was resolved by using double precision 
arithmetic, after it had also been established that the WICS matrices were indeed well-conditioned. 
There is no reason to suspect that such issues will arise again in WICS calculations. Thus the present 
focus is the speed of calculation. 

As a preliminary set of geometrical ideas was developed, the interpretation of iterative solution 
procedure was attempted on small matrices. Preliminary results were not very promising when applied 
to typically large systems however. Computational efficiency was hampered in large matrices due to 
a large number of error-accumulating procedures. The tests finally suggested exploring the current 
linear solver in PMARC [1] in conjunction with the developed geometrical ideas, which resulted in 
improved performance. 
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The questions of solving large linear systems will always arise in our applications since internal flow 
models involving PMARC will typically grid the major portions of a whole wind tunneL However, 
the issues arising from supercomputing and parallel processing will not be attempted here since they 
involve different strategies. 

The basic layout of this report covers both theory and computational experiments. The theoretical 
ideas are substantiated by geometrical interpretations for small matrices. However large matrices 
involving very large solution data files are not presented as outputs. Instead, their performance is 
compared with respect to computational run time using double precision arithmetic. These are 
presented in a tabular fashion and discussed in the results section. The actual fortran programs that 
were developed under this project are available on the sr71 and ra-iris systems at NASA-Ames and 
the author's home directory at Rochester Institute of Technology. Some sample fortran programs are 
attached herewith in the appendix. 


Solution of Linear Systems 

'rhis field of study is very old and practically arises in any area of mathematical interest. The basic 
techniques of linear algebra are learnt in a variety of ways from high school mathematics through 
higher levels involving ideas from topology. The characteristic of all linear systems consists of a 
set of coefficients and unknowns related by a system of equations which never involve any powers 
of the unknowns more than one, neither do any of the equations involve any products of two 
unknown quantities. In a basic appearance we shall call such a set of equations an n-dimensional 
linear system, given by 


a u x 1 + a 12 x 2 +a n x 3 + + a lB x* = b t (1.1) 

a 21 x 1 +a 22 x 2 + a 23 x 3 + + a 2n x n = b 2 (1.2) 

a 31 x, + a 32 x 2 + a 33 x 3 + + a 3n x* = b 3 (1.3) 


a nl x, + a n2 x 2 + a, l3 x 3 + +3™^ = ^ (l.n) 


In the above system, the right hand terms b, through b n are all known, as well as, any terms involving 
the letter "a" with subscripts. The right hand known constants are expressible by a column vector b 
of size n, and the unknown variables x, throughx n may be expressed by another column vector x. Thus 
the linear system may be expressed in a more compact form by A x = b, such that A is a matrix of size 
(n x n) [Le., n rows and n columns]. 

There are various questions of to answer about the existence of solutions and solvability which 
involve the matrix A of different numbers of rows and columns. However in our applications database 
no such cases will arise. Thus we shall always focus on a linear system that poses a unique solution 
and all the solutions arrived at by various means will involve real numbers. Also our system of 
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equations will be non-homogeneous, Le., the vector b will never be a null vector. When a 
homogeneous system of n linear equations are solved, the determinant of A must be zero. Although, 
such systems involving eigenvalues and eigenvectors will be discussed later, the main thrust of the 
solution will involve numerical procedures adopted for non-homogeneous systems. 

The focus of this investigation is in developing and understanding computational approaches. Thus 
traditional solutions of non-homogeneous linear systems given by Cramer’s rule, cofactors and 
evaluation of determinants will be omitted here since these techniques are very expensive 
computationally. The traditional direct solvers include the Gaussian elimination technique, Gauss- 
Jordan technique and the L-U decomposition technique. 


Direct Solvers 


Gaussian Elimination: 

In this approach the idea is to triangularize the matrix A through a set of algebraic reductions such 
that the resulting matrix A' is a strictly upper triangular form. There are different ways to achieve 
this. The most popular technique yields the diagonal elements of A’ as one. The corresponding 
right hand column vector b' after the reductions are utilized in a back substitution process to 
solve for the unknown vector x. The reason Gaussian elimination is claimed as a direct solver is 
because, with sufficient accuracy, the calculations need to be performed only once. With single 
precision arithmetic however, on some limited memory computers, the calculations may need to 
be iterated using an error equation technique. 

Although the process of Gaussian elimination is simple, it does not involve the least amount of 
calculational efforts. It may be shown that the calculations involved in such processes are of the 
order of n 3 operations. There are some associated computational questions to produce robustness 
in such calculations for a general matrix. These involve rearranging the equations so that the 
diagonal elements of A are always the largest. These are called the row pivoting and the column 
pivoting operations. The examples quoted for comparisons later in this report always involved a 
partial row pivoting strategy (for example, see the 9 x 9 truss problem later). 


Gauss- Jordan Elimination; 

This is a step further from the Gaussian elimination process. In this method the eliminations of the 
coefficients are carried out not only below the diagonal but above the diagonal also. The 
advantage of this method is that the reduced matrix A' is a diagonal one, eliminating the need for 
the back substitution process, which is a characteristic of the previous technique. The choice to 
use the Gaussian elimination or the Gauss-Jordan elimination is a personal one since the 
associated procedures involve the same amount of computational efforts. Thus Gauss-Jordan 
technique is not claimable as an improvement over the Gaussian elimination process. For example. 
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the subroutine gss in the current solver in PMARC is a Gauss Jordan one. 


L-U Decompostion: 

The ideas of Gaussian elimination is extended another step in the L-U decomposition process. There 
are different authors for such methods (for example Crout's method, Doolittle’s method, etc).The idea 
is to obtain a simultaneous upper and lower triangularization of the original coefficient matrix A, such 
that A = L.U. There is again no direct advantage of adopting this procedure for computational 
efficiency over the Gaussian elimination technique. However such decomposition ideas are 
worthwhile to understand the direct and iterative numerical techniques presented later. Thus although 
no separate consideration of this technique will be presented in this report, note that the QR 
decomposition technique quoted later to check the ill-conditioning of the PMARC coefficient 
matrices involved such triangularization ideas. 

In summary, the direct solvers available in literature today are all variations of the Gaussian 
elimination process. There are several different methodology to obtain solutions of linear systems by 
direct solvers. However since all direct solvers involve associated computations with the whole 
coefficient matrix A, there is no apparent computational benefit over the original Gaussian 
elimination, since all involve n 3 arithmetic operations. 


Iterative Solvers 

Iterative solvers became more popular with the advent of high speed computers. There are several 
benefits of choosing an iterative solver over a direct solver. Simple iterative solvers invariably 
reduce the programming efforts. However the more sophisticated ones may involve considerable 
amount of complex programming since they can act as a black box for the end user. The 
development of a robust calculation procedure has its roots in the linear algebraic techniques far 
beyond the reach of the end user. We however will- develop some geometrical approaches here in 
an effort to understand the basic iterative techniques. There is a difficulty in presenting 
geometrical ideas beyond three dimensions. The visualization handicap will be supplemented by 
algebraic reasoning. The sections below present the basic and the more recent iterative solution 
procedures. 

To give an example of an iterative process, let us assume an equation: ax 2 + bx + c = 0, where, a, 
b and c are constants and x is the unknown variable. We wish to determine the unknown variable 
x by repetitive calculations called iterations. Note that the equation's solution can be determined 
by the quadratic formula x = { -b ± -/ (b 2 - 4 a c) }/ (2a). However that will be considered a direct 
analytic solution. Instead an iterative procedure may be set up by casting the above equation into 
the following form x(ax+b) = -c, followed by, x = -c/(ax+b). Thus a new value of x may be 
determined from a guessed value of x by the last equation. This is formalized by writing the 
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equation as x (t+1) = -c/(ax a) +b), where (k) represents the iteration number. In this process, the 
guessed solution will be changed and hopefully will be converged to the analytic solution of the 
equation. The theory behind solution of linear equations by iterative processes involves study of 
the procedures and errors associated with iterations and whether or, how quickly the calculated 
solution can be converged to the analytic solution. Note that contrary to this non-linear example, 
our system of equations is linear and associated unknown x is a vector with n components. Given 
below are the popular iterative solvers for linear systems. 


Jacobi Method: 

This technique is the most basic form of relaxation methods and is listed here for reference. The 
actual procedures adopted in this report did not employ this except to calculate part of the solution 
in the successive over-relaxation scheme. The point Jacobi technique involves decomposing A as 
follows. 

The matrix A may be decomposed into a diagonal matrix D, an upper triangular U and a lower 
triangular matrix L. Then the resulting linear system may be written as 


D x^ +1) = -(L + U) x 00 + b , 
where, (k) is the iteration number. 

Thus the solution of the point Jacobi scheme may be sought using 
= . D'l (L + U) x 00 + D ' 1 b 

In the above equation, the matrix, T ; = - D 1 (L + U), is called the associated point Jacobi iteration 
matrix of the linear system A x = b. Note that the second matrix D ' 1 b of the right hand side of the 
above equation is a known static vector, which will not change during the solution process. Whether 
the solution of the point Jacobi method converges to the analytical solution of the linear system 
depends on the norm of the iteration matrix being less than one. It may be shown that this condition 
on the norm is satisfied if the coefficient matrix A is strictly diagonally dominant. 

The speed of calculation in a point Jacobi scheme is slow compared to the processes discussed below 
since the iteration vector x (k+1) gets modified only once in each iteration. 


Gauss-Seidel Method: 

If the same decomposition that was used in the point Jacobi scheme is organized a little differently, 
the speed of calculations can be considerably improved. This method upgrades the vector x 
continuously within a single iteration using the most recent values of the components of x that are 
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available. The resulting scheme is called the point Gauss-Seidel scheme, where, 

x™ = - (D + L) 1 (U) x w + (D + L) 1 b 

where, the iteration matrix is T GS = - (D + L)' 1 (U) . 

The point Gauss-Seidel scheme is convergent if and only if lfT GS ll s 1. This condition may again be 
guaranteed if the matrix A is strictly diagonally dominant. The Stein-Rosenberg theorem relates the 
convergence and the matrix norms of the Jacobi and Gauss-Seidel iteration matrices [2]. 


Successive Over-relaxation Method: 

The best speed of calculations utilizing the ideas of the point Jacobi and Gauss-Seidel schemes may 
be derived by a linear combination of those two solutions. The resulting scheme is called the point 
successive over-relaxation scheme (or, simply S.O.R scheme). Here the linear combination of the 
point Jacobi solution, Xj and the point Gauss-Seidel solution, x<; S is achieved to yield x^r as 

Xsor = toXcs + (1 - Co) Xj 

In the above equation, co is called the relaxation parameter. It ranges in values typically between 0 
and 2. However, higher values are possible to yield a convergent solution. When the value of co is 
less than 1 the process is called under-relaxation and when co is greater than one the process is called 
over-relaxation. Although co is typically a constant, it may be chosen as a variable too. In each 
problem an optimum value of the relaxation parameter may be found analytically as well as by 
performing computational experiments. The iteration matrix of the point S.O.R scheme is given by 

t sor = - (D + coL) 1 [co L + (1 - co) U] . 

Using the substitutions of P = D' 1 L and Q = D' 1 U, we may re-write the above as 
T sor = -(I + coP) 1 [Q + (1 - co) I] 
where, I s the identity matrix of size n x n. 

Two important theorems by Ostrowski and Reich relate the iteration matrices of Gauss-Seidel and 
S.O.R. processes [2]. The important condition which guarantees convergence of the above matrices 
in the complex space is that the component matrices L, U and D be Hermitian and positive definite. 
In the applications of our interest, a real, symmetric and diagonally dominant matrix A would meet 
all the conditions above if the acceleration parameter may be maintained in the range 0 to 2. 
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Conjugate Gradient Method: 


All relaxation processes have the basic idea of reducing the residual vector r [= b - A x] as the 
calculation proceeds from iteration to iteration. There are different procedures defined dependent on 
the search direction to modify the trial vector for x. One of the most successful procedures is called 
the conjugate gradient method. In this process, the trial vector, v is modified as follows: 

v’ = v + t p 

where, t is a scalar, v' is the new trial vector and p is the direction vector in which the new solution 
vector, v’ must be searched. The direction of p must not be chosen perpendicular to the error vector 
because then there will be no improvement in the solution vector v\ If the direction of p is chosen 
same as the error vector the previous iteration step, we get the steepest descent method. If the 
direction of p is chosen as same as a weighted constant factor times the previous step’s error vector 
the resulting procedure is called the simultaneous displacement or the Jacobi iteration procedure. The 
best selection of p comes from satisfying the relationship 

A p 00 . p 0 ^ = p w . A p 01 ' 0 = 0 

This is the basis of the so called conjugate gradient method. In the above relation, the dot indicates 
the inner product of two vectors. The vector p (k) is taken as a linear combination of 0 and p^ 0 . 
In this procedure, the residual vectors and p with Ap form two orthogonal systems, hence the name 
conjugate gradient method. 

The conjugate gradient method is by far the best relaxation method among the traditional iterative 
procedures. The convergence in this process is quadratic and is guaranteed within n iterations where 
n x n is the size of the matrix A. The only drawback is that the matrix to analyze is required to be 
symmetric. There are variations of the conjugate gradient method possible for asymmetric matrices. 
However, as mentioned in the section below and substantiated by the results later, the best scheme 
tested in this investigation is a variant of the Lanczos method. 


Recent Developments: 

As we shall see in the graphical interpretations of the program gsol later, the search of the solution 
vector is facilitated by a choice of an orthonormal basis. QR factorization schemes and procedures 
associated with the Gram Schmidt orthogonalization yield such orthormal bases. However, if the basis 
can be established in an elegant way there are significant computational advantages. 

In any iterative solution the Markov chain converges depending on the eigenvalues of the iteration 
matrix. If the largest eigenvalue is less than one convergence is assured. However the speed of 
convergence is also associated with the closeness of the eigenvalues. These ideas prompted a host 
of procedures [3, 4, 5, 6] related to the determination of eigenvalues. For large matrices the speed 
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of calculation is of prime importance. If there is a way to determine search directions while the matrix 
sizes are relatively small and self-correct the process while progressively lar ger matrices are handled, 
the process will surely be more desirable. The basis of the current solver in PMARC has this structure 
in the subroutine lineq. This procedure is based upon Davidson's method [7] of determination of a 
few of the smallest eigenvalues and eigenvectors of a large matrix. 


Computational Efficiency 

There are several ways computational efficiency is measured in solving linear systems. Most 
computational processes are dictated by the accuracy in calculations. This is very important in 
both direct and iterative solution methods. However accuracy plays different roles in these two 
methods. While direct solvers must be very accurate because calculations are performed only 
once, iterative solvers can handle more errors in early iterations if they can be annihilated rapidly. 
There is also the question of conditioning of matrices. If a matrix is well conditioned, any small 
perturbation in the coefficients will not produce large perturbations in the solution vector. 
Assuming the matrices to explore are well conditioned, errors can only be of one kind - rounding 
errors. Unless any solution scheme models after application of some series, there is no concern 
about truncation errors. Rounding errors can be measured as local and global. Local errors that 
take place during each iteration become of global nature when all iterations cumulatively affect 
calculations. A classic example of sensitivity of a solution scheme on global errors is the Lanczos 
scheme. For quite some time, this elegant method was overlooked because the failures were not 
pegged down to cumulative errors. This important lesson that the author learnt is reflected in the 
design of gsol. 

The speed of calculations is reflected in the number of arithmetic operations that each scheme is 
required to perform. The analytical estimate is typically related to the size of the matrix, n. For 
example, a scheme requiring n 3 operations would take approximately 8 times the computational 
effort if the size of the matrix is doubled. The other measure of speed as discussed before (for an 
iterative process) is going to be reflected in how rapidly the residuals are annihilated. This is 
typically measured by the ratio of the absolute values of the residual norms between two 
successive iteration steps. This can also be assessed by the Euclidean norm of the iteration matrix. 

Another measure of computational efficiency is given by robustness. A robust calculation 
procedure will normally not fail if the conditions of matrix coefficients, the right hand side column 
vector, the starting guess solution, etc. change. This was the focus in the early part of this work. 

Finally, modem linear algebra has been tremendously impacted by the architecture of computers. 
As the demand of speed increases, the concepts are modified to suit the need. Today the product 
of two matrices are done by block concepts much more than the conventional earlier methods. If a 
matrix is symmetric, half data storage is exploited. These and parallel architecture are keys of 
modem computing. This work exploited some storage optimization for 2004 x 2004 matrices 
primarily from the need to save memory on the sr71 machine. Also some modularization was 
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adopted for programming ease. Other concepts of modem architecture were not exploited. 


Graphical Interpretation of Residual Correction Schemes 

In this section we shall explore the geometrical ideas associated with residual correction methods. 
The ideas will be developed from geometrical to algebraic reasoning for higher dimensions. Let A 
x = b be represented for a planar case first. Let a u + a 12 x 2 = b 3 and a n x, + a n x 2 = b 2 represent 
two straight lines in the x-y plane given by OP and OQ, with O as their point of intersection. The 
point O's coordinates are the so called solution of this system of equations which we wish to 
arrive at iteratively. Let the point S represent the coordinates of the starting guess solution. To 
reach the point O from S, let us first drop a perpendicular on OP from S. Let the point of 
intersection be A. Subsequently another perpendicular may be dropped on OQ from A. Let that 
point of intersection be B. With the knowledge of the points A and B we may drop one last 
perpendicular from B on OP to obtain C. Then the solution O of the system may be arrived at 
from A in one movement along AC by the amount AO, where, AO = AB 2 /AC. 



Figure 1. 


This may be shown easily from the similar triangles OAB and ABC where the angle CAB is 
common. Instead of using the direct value from this formula for AO, suppose there is a multiplicating 
factor to used with AO to reach O from A. This may be viewed upon as an accelerating factor in 
higher dimensions. For a set of n equations in n unknowns, 

resl = b, - [a n a 12 a l3 a t J . x$ , represents the residue from the first equation, which also 

represents the perpendicular distance of the point O from the hyperplane: 
a u x, + a 12 x 2 + a 13 x 3 + + a ln x* = b, 

Now this distance can be used to arrive at the point A if a component can be used in the direction 
toward A to obtain the position vector x A . Thus, 


" *5 + 


res] 


l a il *12 *13 *”* *lJ 


7 ^*u *12 *13 — *lj 
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In a similar manner, x B may be written from x A as 


*4 + 


rns2 


Kj ^22 ® Z 3 •— 


2 1*21 °22 °23 a 2jJ 


where, res2 = - [a^ a^ a^J , x A . Similarly Xc may be obtained from x B by projecting into 

the first hyperplane again. After C is obtained a single evaluation of O is made with the relaxation 
factor co. Then the resulting point becomes a new point O and the process keeps repeating. In this 
manner, each new equation is selected from a set of n equations and the projections are made on the 
first hyperplane. This is the basis of the program qmconj given in the appendix. 

An alternate to the above program will be instead of projecting each new residue onto the first 
hyperplane, each two new residues are formed from the sequence of equations (1.1, 1.2), (1.2, 1.3), 
(1.3, 1.4), etc. This process works faster with larger values of <o and is the basis for program simconj 
given in the appendix. A proof of convergence for the above procedures is given below. 

Let Xq* be the new guess point arrived at after one iteration of the qmconj process. Since the vectors 
AB, BC and AC may written as 

ab = x B - x A , abl = V (ab.ab) 
be = x c - x B , acl = V (acac) 
ac = x c - x A , adl = to abl 2 /acl 

Thus, Xo* = x A + a) abl 2 (Xc - x A ) / acl 2 = Q Xc + (1 - Q) x A , where, Q = co abl 2 /acl 2 . 


Note the structure of the acceleration in the last step above emerges exactly like the point successive 
over-relaxation scheme, where the position vectors for points C and A serve the replacement for the 
Gauss-Seidel and Jacobi iterations before. With the matrices substituted for a 4 x 4 system the three 
basic steps to the solution of the new guess vector would be given by 


( a u a \) 


+ M. x 
in 1 ° 


(®2I 


1C 


M 2 X . 


A* 


X C ~ 


( tf ll 


+ M. x. 

in 1 1 


where, the repeated index i indicate summations over i =1,4 and the matrices M, and M 2 are 
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and. 




I.J 1 L.) ^^12 ~«n a x4 

Vu "u®!, V»„ 

~ g H g 12 ^ g 12 ^ ^lAi ~ fl X2 g I4 

* 1*1 **\Pu a ifiu a \Pu 

^nfij (1 fl i3 Vu 

a i/ , u Vu fl u tf u g U®l/ 

- g ll fl 14 • g 12 fl 14 ^1^14 ^ g 14 ^ 

"x^U Vu Vtt V» U 


t °2l V -^2X^2 -*21*23 -*21*24 

fl 2<*2< *21*2 *2<*2» *2**2f 

^1^22 ^ *22 ^ ~*22*23 ~*22*24 

a i/ i Ti *2i*2i *21*21 *21*21 

‘ g *X*23 ~*22*23 ^ *13 ^ ~*23*24 

Va *21*2, V21 ***21 

**2X^24 ~*22*24 ~*23*24 ^ **24 ) 

* U *21 * 21*21 * 21*21 * 21*21 


Similarly, M 3 and M 4 matrices may be written following similar structures as M, and M 2 . Thus the 
final expressions for the three steps Xq*, Xq” and Xq*” before arriving at the upgraded initial guess 
may be derived. Finally the iteration matrix for the program qmconj may be obtained as 
T qm = (1 - Q)“' 1 (M l M 4 M 1 )(M 1 M 3 M,)(M 1 M 2 M 1 ) . In a similar manner the iteration matrix for the 
program simconj may be written as = (1 - G)- 1 (M 3 M 4 M 3 )(M 2 M 3 M 2 )(M 1 M 2 M 1 ). 

Several computational experiments were performed to claim the convergence of the above iterative 
processes. It was shown that in each case the matrices were convergent. The norms of T m and T OT 
were very small indicating rapid reduction of the residuals. Thus the procedures qmconj and simconj 
could be applied for small matrices fairly reliably, each time converging. Another advantage of these 
schemes over other previously discussed iterative methods was they could be applied for ill- 
conditioned matrices (see next section) as well as general solution procedures, where typically Gauss- 
Seidel type methods would fail due to the lack of diagonal dominance and conjugate gradient type 
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methods would fail due to the lack of symmetry. However, the procedures are not very efficient 
computationally specially for large matrices. Thus although the schemes were robust they were also 
computationally expensive. Nonetheless there was an important lesson to be learnt from these 
exercises, which reflected in basic geometrical interpretations of the convergence processes in a 
relaxation method. 

Returning to the planar case of the above methods, we may obtain another arrangement by this 
process to reach O. Instead of obtaining the point A and then using it to obtain point B, let SA and 
SB be two perpendiculars from S onto OP and OQ. Then the arrangement would look like 



Note that the distance SA represents the perpendicular distance of a point from a line. In higher 
dimensions, the corresponding distance represents the perpendicular distance of the point S from the 
n-dimensional hyperplane. If the figure 2 is interpreted another way, one would discover a nice 
geometrical feature imbedded in it. The perpendiculars SA and SO tend to suggest that OS resides 
on the diameter of a circle. In fact the points S, A, B and O are all on this circle with the diameter OS. 
If there was a way to quickly obtain the center C of this circle, obtaining the solution O is just a 
matter of doubling the distance from S to C. This also raises the question that although the above 
geometric reasoning holds for a circle, will it also hold for a sphere in n-dimensions. It turns out that 
the point C would become the circumcenter for the n-dimensional solid. The name circumcenter 
arises from the fact that, for a triangle, the circumcenter is the perpendicular bisector of all the three 
sides, and is at an equal distance from all the three vertices. The process to obtain C in n-dimensions 
was first implemented in the program gsol2v as a direct solver. There are two options to solve for the 
circumcenter. The first one is using the orthonormalization process as mentioned above. An alternate 
procedure is to calculate the Menger's determinant [8]. This process determines the circumradius first. 
However the solution of the resulting system must again be carried out using some speedy iterative 
procedure. Thus the direct solver was later modified as an iterative process in the program gsol to 
obtain the circumcenter using the current solver in PMARC. 

The direct solution of the circumcenter C is found by going to the midpoints of each segment SA, SB, 
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etc. for the n perpendiculars to the hyperplanes for which linear system the solution is sought. This 
resulted in a Gram-Schmidt type orthogonalization procedure. The equations can be cast into a simple 
form starting from the point S. If the base points of perpendiculars to all the hyperplanes are 
obtained, each two new basepoints will make a triangle with the first basepoint of which the 
circumcenter is equidistant from the vertices. However for storage purposes of the large matrices, 
all basepoints are not calculated upfront. First the perpendiculars to the first two hyperplanes, A and 
B are obtained. Starting from A, the midpoint of segment AB is calculated. As each new base point 
C is introduced with the two original basepoints A and B, a scalar multiple X is used to determine the 
proportional distance of the circumcenter on the normal to AB by the formula 

■Hfr/O - (*/*.)] - (* e - XJ'MC 

X - 

(x a -xj-nu 


where, me and nu are the midpoint of segment ab and unit normal to segment ab. The process is 
repeated as new points are read in. The total amount of computational effort is minimized by the 
nested looping in the program. 


Defective Matrices 

As mentioned in the introduction, the difficulties of the PMARC convergence was believed to be 
the ill-conditioning of the coefficient matrices. It was proved later that the WICS calculations will 
always involve diagonally dominant coefficient matrices. Such matrices will not typically have 
small determinants. Thus the small perturbations in the coefficients will not produce large 
perturbations in the solution vector. However, since this research started before the stage when 
the well-conditioning could be claimed, the schemes developed under this research were tested 
with Hilbert matrices. A Hilbert matrix is very ill-conditioned and solutions are almost impossible 
for larger systems. The tests with Hilbert matrices were performed for a maximum of 20 x 20 size. 
It was believed that if a solution scheme successfully performed with a Hilbert matrix, it probably 
was very robust. 

The Hilbert matrix has A's coefficients of the following form: a^ = l/(i+j-l), [Is i,j s n] with the 
right hand vector b varied differently for different computational experiments. Three categories of 
the b vector was tested - i) larger terms toward the top of the column, ii) even size terms in the 
column and iii) larger terms toward the bottom of the column. It was found that the last category 
was the most difficult to solve in most of the cases. Both the routines simeonj and congrad 
performed better than the Gauss Seidel technique when the matrix sizes increased. In this context 
it may also be mentioned that simeonj or qmconj type approach has an additional advantage. 
These procedures will hold even if the linear system is over-determined, meaning that there are 
more number of equations than the number of unknowns. Then although standard reduction 


15 



processes would not work, these procedures will produce a "limit cycle" solution without much 
difficulty. Note that the conjugate gradient method is especially suitable for a symmetric matrix 
even though the matrix is ill-conditioned. In the case of an asymmetric matrix, the conjugate 
gradient method fails. The results are summarized in Table 2 in the appendix. 


Results and Discussion 

As mentioned previously, this work was started as the current solver in PMARC produced 
difficulty in convergence in certain configurations of WICS calculations. The nature of the current 
solver in PMARC was not known to the author. WICS calculations are based upon a procedure 
similar to the error equation approach mentioned earlier in the section discussing direct solvers. 
Thus it was important to test the convergence when the parameter solres was tightened. This 
action typically produced oscillatory behavior and no convergence for certain singularities. The 
first pan of this investigation was therefore to determine the cause of the failure in the WICS 
computations. This was done by extracting a typical coefficient matrix that was utilized during the 
WICS calculations with doublets. This matrix was subsequently subjected to a direct Gaussian 
elimination technique and the QR decomposition technique outside the PMARC set of 
calculations. The prescription of the gridding and section definitions in PMARC were changed 
and a double precision arithmetic was used. From these actions the convergence difficulties were 
overcome. An independent look at the structure of the coefficient matrices showed that the 
diagonal dominance can always be guaranteed in such calculations. The difficulties experienced 
earlier can be ascribed to the loss of orthogonalization typical to a Lanczos process. 

The appendix lists samples of the programs conj, qconj, qmconj, simconj, gsol2v and gsol. The 
results produced by these programs are compared with the Gauss-Seidel iterative technique and 
the conjugate gradient method (see congrad in the appendix) for small matrices. Also tests were 
performed for some sparse (9 x 9) planar truss problems, the ill-conditioned Hilbert matrices and 
some medium size [(79 x 79) and (84 x 84) respectively] external flow calculations over NACA 
0012 and NACA 2412 airfoils. Finally some solutions on the PMARC calculations with a (2004 x 
2004) matrix are presented. Then the comparisons are made of the Gauss-Seidel method, the 
S.O.R method, the conjugate gradient method, the direct solver gsol2v, the direct solver in Gauss- 
Jordan method, the iterative solver gsol and the current PMARC solver. 

First consider the Table 1 in the appendix. The reason 4x4 matrices were chosen for most of the 
testing in this table is because a size 4 x 4 is general enough for the theory to hold even in higher 
dimensions and yet the calculations are not time consuming. Thus many comparisons could be 
made. Also note that the calculations presented in the tables are a subset of the case files available 
in the author's directory with more details. In all calculations in this section, the convergence 
criterion was chosen to be 0.5 x 10 s in the 2-norm of the residual vector. Also smaller matrices 
were run with single precision arithmetic whereas, the large matrix in Table 3 results below were 
all run with a double precision arithmetic. 
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As one can see from Table 1, the number of iterations taken by the examples of 4 x 4 matrices 
simconj and qmconj performed better than the number of iterations taken by the Gauss- Seidel method 
in most of the cases. These routines also performed much better than simple relaxation method conj 
and an early variant of the geometrical approach, qconj. The rapidity with which these calculations 
were achieved could be compared with the conjugate gradient method if proper acceleration 
parameters were selected. Best acceleration parameter choices were dependent on the type of the 
matrices selected. In general, for the 4 x 4 cases the oo selection was around 1.2 - 1.3 for qmconj and 
around 1.6 - 1.7 for the program simconj. Exact values for each case is included in the casefiles. 

The off-diagonally dominant matrices did not yield solutions by the conventional iterative methods 
such as Gauss-Seidel and S.O.R. methods, whereas these two routines consistently produced 
solutions. Also note that qmconj and/or simconj typically take fewer iterations to converge than 
Gauss-Seidel method even for the 4 x 4 diagonally dominant matrices. As the matrix became sparse 
though, the Gauss-Seidel method yielded better performance (see the 9 x 9 planar truss example). An 
interesting case applies to the diagonally dominant Hilbert type 20 x 20 matrix, where the original 
diagonal elements of the Hilbert matrix was modified to maintain diagonal dominance. In this case, 
the best performance was by the use of the conjugate gradient method confirming the superiority of 
this method for symmetric systems. Another interesting case was that of the under- determined system. 
In such situations, the matrix reductions produce a null row suggesting no unique solution of the 
system. Thus conventional methods would not work again. However these geometry based methods 
did not produce a division by zero, and quickly produced a solution satisfying the equations. 

For the external panel method solution applied to the NACA 0012 airfoil, the best performance was 
from using the qmconj matrix among the geometry based methods. For some reason, believed to be 
rounding errors, the program simconj did not perform very well Also note that for the well 
conditioned matrices the Gauss-Seidel method performed better. For the corresponding cases of the 
asymmetric NACA 2412 airfoil, qmconj arrived at the solution with fewer iterations in both 0° and 
90° orientations. In summary, the routines qmconj or simconj achieved the objectives of robust 
calculations typically for the off-diagonally dominant cases when the conventional iterative methods 
had either no solutions or, had difficulty in arriving at them. 

Table 2 quotes the results for some Hilbert matrices. As mentioned before, these matrices were run 
with three different choices of the right hand vector. The reason for this choice was to effect the ill- 
conditioning from mild to severe. In mild ill-conditioning the right hand column vector had larger 
elements toward the top rows and for the severe conditions the larger elements were toward the 
bottom rows. The tight convergence criterion was relaxed a decimal digit (Le., e = 0.5 x lO -4 ) to 
allow faster convergence in the matrix sizes larger than 2x2. Also the larger matrices were run with 
mild ill-conditioning to speed up calculations. The calculations showed an advantage of using the 
qconj and qmconj routines. The routine qconj is clearly the best performer in all computations except 
the case 4 for a 3 x 3 matrix. This success may be attributed to the structure of the cycling in this 
program. All the geometry based routines suffer error accumulation and this routine helps reducing 
it through cycling. As the matrix sizes grew much bigger with ill-conditioning, these routine suffered 
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larger round-off and took too long to converge. This prompted the strategy change and the 
development of the second geometry based solver gsol2v (see Table 3). 

gsol2v is called a direct solver because the determination of circumcenter and the circumradius are 
obtained by direct application of Gram- Schmidt type orthogonalization process. This method was 
sought as a new approach to obtain solution of linear systems. Although the picturization of the 
circumcenter concept is not possible for higher dimensions, the success of the approach shows that 
the low dimensional geometric reasoning can be extended algebraically to higher dimensions. Several 
issues of convexity must be discussed to obtain the algebraic proof of concepts for conventional 
iterative methods. The attempt here was to bypass these by the current geometrical approach. From 
this standpoint, this approach was successful A new correlation between the geometry of spheres and 
solution of linear system was obtained. In planar dimensions, the nine-point circle establishes the link 
between the centroid, the orthocenter and the circumcenter. In higher dimensions, the distances and 
scalar products were needed instead of angles. 

The use of gsol2v required much larger solution time compared to the direct solvers like Gauss- 
Jordan method. It may be mentioned here that among all direct solvers simultaneously involving the 
complete set of coefficients in the matrices, the Gaussian elimination or Gauss-Jordan methods are 
the most economical. Thus it was expected to be slower than those methods. However the argument 
in favor of gsol2v was that it could be applied without pivoting strategies. To improve the speed of 
this solver, the geometrical concepts in this method were combined with an excellent iterative solver. 
This is the idea behind the program gsol. Also note that although the direct solver gsol2v took more 
time than the program impjord , it was less time consuming than both point iterative methods such 
as Gauss-Seidel or the S.O.R. method (with the optimum acceleration parameter * 1.6). 

The algebraic iterative solver used in the lineq routine of PMARC has an excellent structure. It 
requires very little storage and is computationally very efficient. Whereas a traditional Gauss-Seidel 
method took over 15 hours for solving the 2004 x 2004 matrix, this solver obtained the solution in 
4 minutes. The secret of this solver drew the author to study the eigenvalue based methods after 
Nesbet, Shavitt, Bender and Davidson [3-6]. There are a set of programs nesbet, nespy, nesl, nes2, 
etc. and the resulting output files available for these methods in the authors directories. These 
programs calculate the lowest eigenvalues and eigenvectors for large symmetric and asymmetric 
matrices. However the current solver fills the voids suffered in all those earlier programs. All of those 
[3 - 6] started the search for eigenvalues from the Raleigh quotient approach. However considerable 
difficulty was experienced trying to simultaneously change all components of the search vector. The 
current program can also handle asymmetric matrices. It starts from small matrices and applies 
corrections to the search direction as the matrix grows. It is also fully optimized for large matrix 
calculations. Note that although it uses the Gauss-Jordan method internally to invert matrices, the 
calculations take practically no time (see the impjord results in Table 3) since the main time 
consuming calculations are performed on much smaller matrices. The design of this procedure even 
offers advantages in restarting the unconverged solutions. Combining the geometrical approach in 
gsol2v therefore with this solver yielded the apparent benefits. The resulting process took just one 
more iteration to converge the 2004 x 2004 matrix. Also the calculation time to converge was 
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reduced from more than 3 hours for gsol2v down to 4 minutes in gsol. 

Note the contrasted results by the conjugate gradient method in Table 3. Since the PMARC matrices 
are not symmetric for WICS calculations, conjugate gradient method could not be applied directly. 
The program mi.xcon therefore modified the solution approach by premultiplying the system A x = 
b by A t . This resulted in a symmetric system which could then be solved by the conjugate gradient 
method. This approach yielded the results of the matrix in 1 hour and 47 minutes and using 99 
iterations. However if we subtract the time required to obtain the products A T A and A T b, the basic 
conjugate gradient process took only 12 minutes. Thus calculations of such large matrices proved to 
be extremely sensitive to arithmetic operational counts. This also shows that besides the current 
solver, the next best approach in solving large systems will perhaps have to be based upon the 
conjugate gradient method. The program gsol offers the quickest alternate solution (in a single step) 
as long as the circumeenter could be obtained quickly. Further expansions are possible for this idea. 


Conc lu di n g Re ma rks 

The current investigation was prompted by a difficulty in convergence of PMARC matrices when 
applied for WICS calculations. The objective was to develop a robust solver that would work 
typically when other iterative methods failed to produce solutions. This objective was fulfilled by 
developing a geometry based solver in the approaches of conj, qconj, qmconj and simconj. There 
was an alternate geometrical approach developed parallel to the conventional algebraic 
approaches for iterative methods in the programs gsol2v and gsol. Throughout this investigation, 
an attempt was made to re-visit geometrical topics as a means to enhance physical and intuitive 
understanding of the convergence processes. The most recent and more complex methods are 
perhaps more sophisticated computationally. However a renewed exploration of geometrical 
concepts probably contributes to a much better understanding than abstract ideas offered in them. 
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Appendix 

This appendix contains tables used for comparing results and samples of 7 fortran programs conj, 
qconj, qmconj, simconj, gsol2v, congrad and gsol. The solvers for Gaussian elimination, all of the 
eigenvalue based methods, lineq, QR Decomposition method, Gauss-Seidel and S.O.R. methods 
are not presented here for compactness and are available in the author's directories. Also all the 
details of the casefiles reported in the tables below may be found there. 
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Table 1 

Comparison of Various Methods for Small Matrices 

In the following results, the matrices are specified by their coefficients. The strategy was to 
compare the number of iterations to convergence and typically choose off-diagonally dominant 
matrices so that the robustness of calculations could be tested. Most cases started with the same 
starting solution: x, = x 2 =....= x„ = 1.0 (unless specifically mentioned in the casefiles). The 
convergence criterion was 0.5 x 10' 5 in the 2-norm of the residual. Also the acceleration 
parameters used in these were near optimal and mentioned in the casefiles. 



Coefficients of A and b 
(Comments) 

No. of Iterations to Converge Method 
conj qmconj simconj G. S. C.G. 


2, 1 ,-3, 1 ; 1 ,2,5,- 1 1 , 1 , 1 ,4;2,-3,2,-5 
l,7,5,-4 (off-diagonal dominance) 

n.a. 

28 

12 

failed 

failed 

D 

l,l,l,l;2,l,-l,l;-l,2,3,l;3,2,-2,-l 
2,2, 1,8 (similar to above) 

208 

38 • 

70 

failed 

n.a. 


-1,-4, 2,l;2,-l,7,9;-l,l,3,l;l,-2, 1,-4 
-32,14,11,-4 (mild off-diag dom.) 

42 

15 

9 

failed 

failed 

■ 

l,l,l,0;-3,-17,l,2;4,0,8,-5;0,-5,- 

2,1 

3, 1,1,1 (mild off-diag. dominance) 

97 

27 

33 

failed 

failed 

n 

2,-l,0,0;-l,3,-l,0;0,-l,3,-l;0,0,-l,2 
-1, 4,7,0 (diagonally dominant) 

30 

13 

6 

15 

5 

D 

2,- 1 ,0,0;- 1 ,4,- 1 ,0;0,- 1 ,4,- 1 ;0,0,- 1,2 
3,5, -15,7 (diagonally dominant) 

30 

13 

6 

10 

5 


2,-2,- 1 ,3 ;- 1 ,0,- 1 , 1 ;3, -6,2,4; 1 ,-4,3, 1 
5, 0,7,2 (under-determined system) 

22 

8 

5 

failed 

failed 


l,l,l,l;l,l,-l,l;-l,2,3,l;3,2,-2,-l 
3,2, -2,- 1 (off-diagonal dominance) 

52 

24 

33 

failed 

failed 

9x9 

Truss Problem (no re- 
arrangement) 

148 

117 

145 

failed 

n.a. 

9x9 

Truss Problem (eqns. re-arranged) 

139 

127 

53 

16 

n.a. 

■ 

Re-arranged Hilbert Matrix 

n.a. 

large 

large 

n.a. 

2! 
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79 x 
79 

NACA 0012 Airfoil (Panel 
Method for External Flows) 

568 

72 

400 

44 

n.c. 

79 x 
79 

NACA 0012 Airfoil (Panel 
Method with a =90°) 

906 

83 

938 

53 

n.c. 

84 x 
84 

NACA 2412 Airfoil (Panel 
Method) 

856 

103 

880 

n.a. 

n.c. 


n.a. = not available, n.c. = non-convergent 


Table 2 

Comparison of the Methods applied to Hilbert Matrices of different sizes 


Size 

(Comments) 

Iterations to converge method to specified criterion 
G.S. conj qconj qmconj simconj 

2x2 type 1 

37 

624 

13 

5 

18 

2x2 type 2 

34 

559 

17 

7 

16 

2x2 type 3 

36 

697 

21 

7 

20 

wem 

199 

n.c. 

> 10,000 

49 

104 

■ 

i 

126 

n.c. 

14 

203 

109 

6x6 

427 

n.c. 

27 

129 

1048 

10 x 10 

261 

n.c. 

70 

75 

423 

20x20 

515 

n.c. 

52 

512 

1722 


n. c. means no convergence was achieved in 10,000 iterations. 
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Table 3 

Comparison of Various Methods for a PMARC Matrix of Size 2004 x 2004 


Program Name 

Method Used 

Number of Iterations 
to convergence 

Solution Time 

gssd 

Gauss-Seidel Method 

6936 

15 hours 1 1 min. 

sor 

S.O.R. (co = 1.6) 

5554 

6 hours 22 min. 

gsol2v 

Finding Circumcenter 

Direct Solver 

3 hours 37 min. 

mixcon 

Conjugate Gradient 

99 

1 hour 47 min (*) 

impjord 

Gauss-Jordan Meth. 

Direct Solver 

1 hour 55 min. 

gsol 

Circumcenter Method 

41 

~ 4 minutes 

lineq 

Davidson’s Approach 

40 

~ 4 minutes 


* The time taken by the iteration portion only was 12 minutes. 


Program Listings: 

1) coni.f 


program cooj 
parameter n = 20 

dimension a(n,n), b(n), x(n), y(n\ res(n) 
OPEN(UNIT= 1 3 3LE= , a.dat^STATUS= , OLD , ) 
OPEN(UNIT=14J 1 * * * * * 7 ILE= , b.dat',STATUS= , OLD') 
read (13,*) ((a(ij), j = l,n), i = l,n) 
read (14,*) (b(i) t i = I,n) 
c write<6,3)((a(ij), j=l,n),i=l,n) 

3 formatCU 9f8.3) 
write<6,3) (b(i), i=l f n) 

c 

c The above line format must be changed 

c 

write(6,10) 

10 format(// Results using the Program conji //) 

c initial guess 

itr=0 

do 1 i=l,n 
c 1 x(i)=0.0 

1 x(i)=l.O 

101 k=l 

100 continue 

sum=0. 
ysq=0. 
do 4 i=l ( n 

4 y(i)=a(k i i) 
do 5 i = l,n 
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ysq=ysq+y(i)*y(i) 

5 sun*=sunH-y(i)*x(i) 
res(k)=bOc)-sum 
do 6 i=l,n 

6 x(i)=x(i)-fres(k)*y(i)/ysq 
if(k.eq.n) goto 8 
k=k+l 

go to 100 

c 8 if(itr.eq.l000) wrile{6,*) (x(i),i=l,n) 

c 8 write(6,*) (res(i),i=l,n) 

c write(6 1 12) (x(i),i=l,n) 

8 resmx=0. 

do 11 i=l, n 

resmx=resmx+res(i)*res(i) 

1 1 continue 
resmx=sqrt(resmx) 
itr=itr+l 

if(itr.lLlO) write{6,*) itr, resmx 

if(resmx.gt 0.00005 and. itr.ltlOOOO) go to 101 

write(6,I2) (x(i),i=Un) 

12 fotmat(2 ^'solution: \5fl2.5) 
write(6,13) itr 

1 3 fonnatC The above was obtained in ',i6/ iterations') 
stop 

end 

2) qconi.f 

program qconj 
parameter n= 20 

dimension x0(n), xa(n),xb<n),xc(n),ab(n),bc(n),ac<n),a(n4i) l b(n), 

1 y(n) > res(n) t err(n) 

OPEN(UNIT=13, FrLE^a.dat'.STATUS^OLD') 
OPEN(UNIT=14, FHJE=’b.dat,STATUS=OLD') 
read(13,*)((a(ij)j=l,n)4=l,n) 
read(14 t *)(b(i),i=l,n) 
c write(6,3)((a(ij)j= 1 ,n),i= 1 ,n) 

3 format(lx,9f8.3) 

c write(6,3)(b(i),i=i,n) 

c 

write(6,21) 

2 1 format(//5 x/Results with Program qconj T//) 

c 

c a and b store the original Ax = B 

cc initial guess 

om=1.6 

c om=1.333 

c om=1.6 

itr=0 

do 1 i=l,n 

1 x0(i)=0.0 

c 1 x0(i)=1.0 

c 

c if n is even, proceed exactly. If n is odd, set last cycle steep descent 

c 

1st = n/2 
iodd=n - 2*ist 
write(6,2) ist,iodd 

2 fonnatC ist, iodd=’,2i5) 

100 continue 

do 200 ic=l,ist 

k=2*(ic-l)+l 

kpl=k+l 

c write(6,*) itr, k, kpl 

do4 i=l,n 

4 y(i)=a(k,i) 
sum=0. 
ysq=0. 
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do 5 i=l,n 
ysq=ysq+y(i)*y(i) 

5 sum=sum+y(i)*xO(i) 
res(k)=b(fc>-sum 

c wnte(6,*)(x0(i),i=l ( n) 

do 6 i= 1 ,n 

6 ju(i)=xO(i)+res(k)*y(i)/ysq 
do 41 i=l,n 

41 y(i)=a(kpl,i) 
sum=0. 
ysq=0. 

do 51 i=l,n 
ysq=ysq+y(i)*y(i) 

5 1 sum=sunH-y(i)*xa(i) 
res(kp 1 )=b(kp 1 Vsum 

c write(6,*)(xa(i),i==l,n) 

do 61 i=l,n 

6 1 xb(i)=xa(i)+res(kp 1 )*y(i)/ysq 
do 42 i=l,n 

42 y(i)=a(k,i) 
sum=0. 
ysq=0. 

do 52 i=l,n 
ysq=ysq+y(i)*y(i) 

52 sum=surofy(i)*xb(i) 
res(k)=b(k)-sum 

do 62 i=l,n 

62 xc(i)=xb(i)+resOO*y(i)/ysq 

c calculate ab, be, and ac (vectors & lengths) 

do 43 i=l,n 
ab(i)=xb(i)-xa(i) 
bc(i)=xcO>xb(i) 

43 ac(i)=xc(i)-xa(i) 
sum=0. 
ysq=0. 

do 53 i=l,n 
ysq=ysq+ac(i)*ac(i) 

53 sun«uim-ab(i)*ab<i) 
abl=sqrt(sum) 
acl=sqrt(ysq) 
adl=om*abl*abl/aci 
do 63 i=Kn 

63 xO(i)=xa(i)+adl*ac(i)/ad 

c 46 do 47 i=l t n 

c 47 xO(i)=xb(i) 

c 

200 continue 

c 

c now store the odd-balls 

c 

if(iodd.eq.l) then 
do 44 i=l,n 

44 y(i)=a(n,i) 
sum=0. 
ysq=0. 

do 54 i=l,n 
ysq=ysq+y(i)*y(i) 

54 suiw=sun>+ y(i) *x0(i) 
res(n)=b(ii)-sum 

do 64 i=l,n 

64 xO(i)=xO(i)+res(n)*y(i)/ysq 
end if 

c 

c Error Calculation 

c 

do 65 i=l»n 
sumaO. 
do 55 j=l,n 
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55 sum?=surrn-a(ij)*xO(j) 

65 err(i)=abs(b<i)-sum) 

c write(6,56)(err(i),i=l,n) 

56 formate Eirors:',5fl 0.5) 
emoc=0. 

do 66 i=l,n 

66 entot=errtoC+err(i) *err( i) 
errtot=sqrt(errtot) 

if(entotit.0. 00005 .or. itr.gt 100000) go to 67 
c write(6,*) itr, errtot 

itr=itr+l 
goto 100 

67 write(6,68)(x0(i),i=l,n) 

68 formate Solution:', 5f 11. 5) 
write(6,69) itr.om 

69 format (5 x,’Obtained in \i5/ iterations, with om= f ,fl0.3) 
stop 

end 


3) qmcQniJ 

program qmconj 
parameter n= 20 

dimension x0(n), xa(n),xb(n),xc(n),ab(n),bc(n) ) ac(n),a(iMi),b<n), 

1 y(n),res(n),err(n) 

c real *8 sum,ysq,om,abl > ac!,adl 

OPEN(UNIT=13, FTLE= , a.dat\STATUS= , OLD) 

OPENCUNIT=:l4, FILE= , b.dat',STATUS=OLD') 
c OPEN(UNIT= 15. Fn^xin.dat , JK>RM='FORMATTED , ,STATUS='OLD0 

read(13,*)((a(i j) 1 ,n)4=l,n) 
read( 1 4 , *)(b(i),i= 1 ,n) 
write(6,3)((a(ij)j= l ,n),i= 1 ,n) 

3 fonnat(lx,9f8.3) 
write(6,3)(b{i),i=l ,n) 

c 

write<6 t 2) 

2 format(//5x,’ResuJts using the Program qmconj fif) 

c 

c a and b store the original Ax * B 

cc initial guess 

case4 om=1.22 

chilb4 om=1.97 

easel om=1.68 

ccqm om=1.29 

on^l.3 
epan om=1.75 

c ont=1.0 

itr=l 

c read( 1 5, *)(x0(i),i= 1 ,n) 

do 1 i=l,n 
1 x0(i)=0. 

c 1 x0(i)=l. 

c write<6,3)(x0(i)4= 1 ,n) 

c 

100 continue 

do 200 k=2,n 
c write(6,*) itr, k 

do4i=l,n 

4 y(i)=a<l,i) 
sum=0. 
ysq=0. 

do 5 i=l,n 
ysq=ysq+y(i)*y(i) 

5 sum=sunH»y(i)*x0(i) 
res(l)=b(l>-sum 

c write(6,*)(x0(i),i=l,n) 

do 6 i= 1 ,n 
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c 


xa(i)=xO(i)+res( I )*y(i)/ysq 


6 

do41 i=U 

41 y(i)=a(k,i) 
sum=0. 
ysq=0. 

do 5 1 i=l,n 
ysq=ysq+y(i)*y(i) 

5 1 sunt=sunH y(i)*xa(i) 
res(k)=b(k)-sum 

c write(6,*)(xa(i),i= 1 ,n) 

do 6 1 i=l,n 

61 xb(i)=xa(i)+res(k)*y(i)/ysq 

c write(6 t *)(xb{i),i= I ,n) 

do 42 i=l,n 

42 y(i)=a(U) 
sum=0. 
ysq=0. 

do 52 i=l ,n 
ysq=ysq+y(i)*y(i) 

52 sum=sum+y(i)*xb(i) 
res(l)=b(l)-sum 

do 62 i= 1 ,n 

62 xc(i)=xb(i)+res(l)*y(i)/ysq 

c write(6, *)(xc(i),i= 1 ,n) 

c calculate ab, be, and ac (vectors & lengths) 

do 43 i=l,n 
ab(i)=xb(i)-xa(i) 
bc(i)=xc(i)-xb(i) 

43 ac(i)=xc(i)-xa(i) 
sum=0. 
ysq=0. 

do 53 i=l,n 
ysq=ysq+ac(i)*ac(i) 

53 sum=sunvt-ab(i)*ab(i) 
abl=sqrt(sum) 
acl^sqrt(ysq) 
adl=om*abl*abl/acl 
do63 i=l,n 

63 xO(i)=xa(i)+adl*ac(i)/ad 

c write(6,*)(x0(i),i=l,n) 

c 

200 continue 

c 

c Error Calculation 

c 

do 65 i=l,n 

sum=0. 

do55j=l,n 

55 sum=sum+a(ij)*x0(j) 

65 err(i)=abs(b<i)-sum) 

c write(6,*)(x0(i),i=l,n) 

c write(6,56)(err(i) r i=l,n) 

56 formate Errors:\5f 10.5) 
errtot=0. 

do 66 i=l,n 

66 errtot=errtot+err(i)^err(i) 
errtot=sqrt(emot) 
if(itr.lLlO) write(6,*) itr, errtot 
if(entot.it.0.00005 .or. itr.gL 100000) go to 67 

c if (errtot !t 0.000005 .or. itr.gt 1 00000) go to 67 

if(errtotgt 10000.) go to 67 
itr=itr+ 1 

if(mod(itr,1000).eq.0) write(6,*) itr, errtot 
go to 100 

67 write(6,68)(x0(i),i=l,n) 

68 formatC Solution:', 5 fl 4.5) 
wnte(6,69) errto^i^om 
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69 format(3x, r Err= , J > 10.6,’ Obtained iterations, ont='J8.3) 
stop 
end 


4) simcmf 

program simconj 
parameter n= 20 

dimension x0(n), xa(n),xb(n),xc(n),ab(n) 1 bc(n),ac(n) l a(n,n),b<n), 

1 y(n)jes(n),err(n) 

c real *8 sum,ysq,abl,acl,adl,bcl 

OPEN(UNIT= 1 3, FILE= , a.dat’,STATUS= , OLD') 

OPEN(UNrr= 14 , FILE= , b.dat' > STATUS=’OLD , ) 
c OPEN(UNIT= 1 5, FILE='xin.dat r K)RM=TORMATTED',STATUS='OLD') 

read( 1 3, *)((a(i j) j= 1 ,n),i= 1 ,n) 
read(14,*)(b(i),i=l t n) 
c write(6 ( 3)((a(i j) j=sl ,n),i=l ,n) 

3 fonrat(lx,9f8.3) 

c write(6 > 3)(b<i),i=l ) n) 

c 

write(6,2) 

2 format(//5x,’Results using the Program simconj.f//) 

c 

c a and b store the original Ax = B 

cc initial guess 

ccgm om=1.2 

om=1.6 
cbl2 oro=1.6 

cpan om=.05 

c on^l.2 

c5 on>=1.3 

c on*=l.78 

chilb4 ocip=1.95 

chilb3 on^l.89 

cas©4 om=2.9 

c20tr ont=2.368 

do 1 t=t,n 
c 1 x0(i)=0. 

1 xO(i)=l, 

c read( 1 5,*)(x0(i),i= 1 ,n) 

itr=l 
c 

100 continue 

do 200 k=l,n-l 
kpl=k+l 

c write<6 t *) itr, k 

do 4 i= 1 ,n 

4 y(i)=a(k»i) 
sum=0. 
ysq=0. 

do 5 i=l t n 
ysq=ysq+y(i)*y(i) 

5 sum=surm-y(i)*t0(i) 
res(k)=h(k)-sum 

c write(6,*)(x0(i),i=l,n) 

do 6 i= l ,n 

6 xa(i)=x0(i)+res(k)*y(i)/ysq 
c 

do 4 1 i= 1 ,n 

41 y(i)=a(kpl,i) 

sum=0. 
ysq=0. 
do 5 1 i= 1 ,n 
ysq=ysq+y(i)*y(i) 

5 1 sum=sun>+y(i)*xa(i) 

res(kp 1 )=b(kp 1 )-sum 
c write(6,*)(xa(i),i=l,n) 
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do 61 i=t,n 

61 xb(i)=xa(i)+res(kpl)*y(i)/ysq 
do 42 i=l,n 

42 y(i)=a(k,i) 
sum=0. 
ysq=0. 

do 52 i=l,n 
ysq=ysq+y(i)*y(i) 

52 sum=suim-y(i)*xb(i) 
resOc)=b(k)-sum 

do 62 i=l,a 

62 xc(i)=xb(i)+res(10*y(i)/ysq 

c calculate ab, be, and ac (vectors & lengths) 

do 43 i=l,n 
ab<i)=xb<i)-xa(i) 
bc(i)=xc(i)-xb<i) 

43 ac(i)=xc(i)-xa(i) 
sum=0. 
ys<f=0. 

do 53 i=l,n 
ysq=ysq+ac(i)*ac(i) 

53 sum=sum+ab(i)*ab<i) 
abl=sqrt(sum) 
ac!=sqrt(ysq) 
adl=om*abl*abl/acl 
do 63 i=l,n 

63 xO(i)=xa(i)+adl*ac<i)/acl 
c 

200 coatinue 

c 

c Error Calculation 

c 

do 65 i=l,n 
sum=0. 
do 55 j=l,n 

55 sum=suro+a(ij)*xO(j) 

65 err(i)=abs(b(i)-sum) 

c write(6,56)(eiT(i),i=l,n) 

56 formate Errors:', 5fl 0.5) 
errtot=0. 

do 66 i=l,n 

66 errtot=errtot+err(i)*err(i) 
errtot=sqrt(errtot) 

if(errtotlL0.00005 .or. itr.gt 1000000) go to 67 
if(errtoLgt 10000.) go to 67 
if(itr.ltlO) write(6,*) itr, errtot 
if(mod(itr,1000).eq,0) write(6,*) itr,entot 
itr=itr+l 
go to 100 

67 write(6,68)(x0(i)4=l,n) 

68 formatC Solution:’, 5fl 1 .5) 
write(6,69) errtot,itr,om 

69 format(5x,’Err='J10.6,’ Obtained in \i7,’ iterations, with otrt=’/10.3) 
stop 

end 

5) conprad.f 

program congrad 
parameter n = 20 

dimension a(n,n), b(n), e(n), r0(n),pl(n), pk(n), rl(n),x(n) 
dimension ap(n) 

open (unit = 3, file = 'a.dat\ status = *ol<f) 
open (unit = 4, file = ’b.dat', status = ’ol<f) 
read(3,*)((a(i j), j= 1 ,n),i= 1 ,n) 
read(4,*)(b(i),i= 1 ,n) 
write(6,ll 1) 

1 1 1 format^/results obtained by congradf: /) 
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do 100 it= 1,1 0000 
if(itge.2) the a 
sumrO = 0. 
sumrl = 0. 
do 12 i = l,n 

sumrO = sumrO + rO{i)*rO(0 
sumrl = sumrl + rl(i)*rl(i) 
ekml = sumrl /sumrO 
do 13 i = l, a 

pk(i) = -r l(i)+ekml *pk(i) 
endif 

write(6**) ekml 
wnte(6,*) (pk(i),i=l,n) 

do 2 i=l,n 
sum = 0. 
do 3 j=l,n 

sum = sunn-a(i j)*x© 
r0(i) = sum + b(i) 
if(iteq.l) pk(i) = - rO<i) 
pl(i) = - rO(i) 
do 4 i=l,n 
sum = 0. 
do 5 j =l,n 

sum = sum + a(ij)*pk(j) 
ap(i) = sum 
write(6,*) (ap0),i=l,n) 
sumr = 0. 
sump = 0. 
do 6 i =l,n 

sumr = sumr + rO(i)*TO(i) 
sump = sump + ap(i)*pk(i) 
write(6,*) sump 
qk = sumr/sump 
do 7 i = l, a 
x(i) = x(i) + qk*pk(i) 
rl(i) = r0<i) + qk*ap(!) 
write(6,*) qk 
write(6 ( *) (x(i),i=l,n) 
write(6,*) (rl(i)a=l,n) 
sumr = 0. 
do 8 i=l,n 

sumr = sumr + rl(i)*rl(i) 
write(6,*) sumr 
wrile(6,108) 

formatC 7 

sumr = sqrt(sumr) 
itr = it 

if(mod(itr,1000).eq.0) write(6,*) itr^umr 
if(itr.le.l0) write(6,*) itrsumr 
write^,*) itr .sumr 
if(sumr.lt0.000005) go to 101 
continue 
go to 103 
write(6,102) itr 

focmat(3 x, 'solution converged in\i4,’ iterations') 
write(6,105) (x(i),i=l,n) 
go to 106 
write(6,104) itr 

format(3x,'Solution did not converge in',i6,' iterations') 

format(5fl4.6) 

stop 



6) £sol2vf 


c 

c 


111 

l 

c 

c 


3 


c 

c 

c 


2 


c 

c 

c 

c 

c 


8 


10 

c 

c 


19 


c 


11 

c 

c 


12 


13 


c 

c 

c 


program gsol2v 

parameter n = 2004 

dimension b{n), u(n),v(n,n),p(n),po(n) 

dimension x(n),a(n),tn(n)il (n)42(n),ab(n),bc(n) 

x(i), p(i) are the 1 plus n base projection points 

open (unit = 3, file = 'a.dat', status = ’old*) 
open (unit = 4, file = Udat', status = ’olcf) 
open (unit = 7, file = ’gsol2v.out', status = 'qcw') 
read(4,*)(b(i),i=I,n) 
write(7,l 11) 

format^results obtained by gsd2v.f:7) 
do 1 i =l,n 

x(i) = lVsqrt(float(n+l)) 
p(0 = *(i) 

do 9 i=l,n 

read(3,*)(a(j) j= 1 ,n) 
do 3 jsl^i 
C(j)=xO) 

call doCpta^iLsun) 
call dotp(&,a,ousu) 
rs=b(i>sun 

aeate p(n) to store the n base points, ooe at a time 

do2j=U 
po(j)=P(j) 
do7 j=l,n 

p(j)=rs*a(j)/su + x(j) 
write(6,*)(p(j) j= 1 ,n) 

Start the nested Looping 

Use the variable v(n,n) to store the orthooormal base 

if(i.eq.l) then 

do 8 j=l,n 

*b(j)=p(j)-x(j) 

call dotp(ab,ab,iuuni) 

abl=sqrt(sum) 

rl=abl/2. 

do 10j=l,n 

u(j)=ab(j)/abl 

tn(j)=x(j)+rl*u© 

calculate al and use later ! 
do 19 j=l,n 
fl(j)=x© 

G<j>xO) 

call dotp(fl42,n,al) 
doll j=l,n 

v(ij)=u(j) 

else 

j greater than 2 below 
do 12 j=l f n 
bc(j)=p(j>-poO) 
call docp(bc,bc,n^um) 
do 13 j=l,n 
u(j)=bc(j)/sqrt(sum) 

write(6,*)' u =\(u(j)j=Ln) 

Gram-Schmidt Orthonormalization 



n n 


c 


do 14 k=l,i-l 


c 

15 

16 
14 


17 

18 
c 

c 

c 

c 

20 


21 

c 


22 


23 

c 

9 

c 

201 

c 

c 

207 


1 


7) fsol.f 


do I5j=l,n 

fl(j)=u(j) 

f2(j)=v(kj) 

call dotptfl j7,n,sum) 

do 16 j=l,n 

v(ij)=v(ij)+sum*v(kj) 
cooiiaue 
do 17 j=l,n 
v(ij)=u(j>v(ij) 

fl(j)=v(ij) 

f20)=v(ij) 

call dotp(fl 42,n^um) 

do 18 j=l,n 

v(ij)=v(ij)/sqrt(sum) 

write(6,*)' v= ,(v(ij)j=l ( n) 
Proceed with the new vector 

do 20 j=l,n 

fi(j)=p(i) 

f2(j)=p(j) 

call dotp(fU2.n,bei) 
do 21 j=l,n 
n(j)=PO>-xO) 

f2(j)=tn(j) 

write^,*)' tns',(tn(j)jal ( o) 

call dotp(f 1 42,n,gam) 

rnu=(bet-al)/2, - gam 

do 22 jal.ii 

G©=v(ij) 

call dotp(fljf2,iude) 

rl=rnuAde 

do 23 j=l,n 

tn(j)=tn(j)+rl*v(ij) 

write(6,*y tn=',(tn(j)j=l>n) 

end if 

continue 

do 201 i=l,n 
x(i)=x(i)+2.*(m(i)-x(i)) 

write(7,*) (x(i)>i=l,n) 
write(6,207) (x(i),i=l,n) 
format (2x,5fl 5.5) 
stop 
end 

subroutine doCp(a,b,n,sum) 

dimension a(n),b(n) 

sum=0. 

do l i=i, n 

sum=sunn-a(i)*b(i) 

return 

end 


Program gsol 
Parameter n = 2004 
Include param-dal' 

Include commotU’ 

Dimension x(n),a(n) J2(n),p(n) 

X(i), p(i) are the 1 plus n base projection points 


Open (unit = 3, file = a.daf, status = 'old) 
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Open (unit = 4, file = b.dat’ , status = ’old*) 

Read(3, *)((dubicww(i j) j= [ ,n),i= I ,n) 

Read(4, *)(rhsv(i),i= 1 ,n) 

Write(6,l 1 1) 

1 1 1 Format(//results obtained by gsd22.f:7) 

Do 1 I =l,n 
i X(i) = 0. 

C 

Do 9 i=l,n 
C 

C Read(3,*)(a(j) j= 1 ,n) 

Do 12 j=l,n 

12 A(j)=dubicww(iJ) 

Call dotp(a,x,n^un) 

Call dotp(a,a,n,su) 

Rs=rhsv(i)-sun 

C 

C Create p(n) to store the n base points, one at a time 

C 

Do 7 j=l,n 

7 P(j)=rs*a(j)/su 
Do 8 j=l,n 

8 F2(j)=p(j)+x(j) 

Call dotp(f2/2,n,sup) 

Call dotp(x,x,n,sux) 

C Write(7, *)(p0) 1 ,n),(sup-sux)/2. 

Rhsv(i)=(sup-sux)/2. 

Do 9 j=l,n 
Dubicww(iJ)=p(j) 

9 Continue 
C 

Call gsd22 
C 

Do 10 i=l,n 

10 X(i)=2.0*dub(i)-x(i) 

Write(78,l 1) 

1 1 Focmat(5x,'with the above circumcenter, the calculated solution is: 1 ) 
Write(78,*) (x(i),i=l,n) 

Stop 

End 

C 

Subroutine dotp(a,b,n,sum) 

Dimension a(n),b(n) 

Sum=0. 

Do 1 i=l,n 

1 Sum=suro+a(i)*b(i) 

Return 

End 

♦deck doublet 

Subroutine gso!22 
C 

C- 

C 

Include param.dat' 

Include 'common/ 

C 

Open(unit=16, rde='data6'ioni^' formatted 1 ,status=’aew') 
Open(umt=20 1 forn^bnformatted , ,status='unknown') 

C 

Do 2 I = Cnspdim 

2 Diag(i) = dubicww(i,i) 

C 

C Rewind all scratch file 20 and assign unit number 
C 

I mu = 20 
Rewind imu 
C 
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C Check to make sure that if inram is not 1, that it is set equal to nspdim 
C 

If(i mam. ne.l) then 
If(i mam. ne. nspdim) then 
Write( 16,601) 

Stop 

Endif 

Endif 

C 

C Start the time step loop 
C 

Tstime = 0.0 
C 

C Write input data to output file 
C 

Write( 16,603) 

C 

603 format( l x, 1 0(/),3 1 x,57C*)// 

+ 54x, program pmarcY 

+ 4 5x, 'release version 12.21: 03/04/947/ 

+ 5 1 x,’matrix solver extracted from pmarc’/) 

C 

Call solver 
C 

Open(umt=78, file=’gsol.out\ status ='nev/) 

Write<78,*) (dub(i),i= l, nspdim) 

C 

C Close and delete the scratch files 
C 

Close(u nit=20,status=’delete') 

C 

Return 

600 fcrmat(//lx, , time step\i4) 

601 format(//l x, parameter imam not set to l or nspdim in pmarc/ 

+ 1 x, 'source code. Reset this parameter, recompile code’/ 

+ 1 x/and try again. 1 ) 

End 

C 

♦deck solver 

Subroutine solver 
C 

C 

c 

Include paramdat’ 

Include 'common.f 
C 

C Update the pdub array so that it always holds the previous step’s doublet 
C Solution 
C 

ltstep=0 
Npan — • nspdim 
C 

Do 10 i=l,npan 
If(itstep.eq.O)then 
Pdub(i) = rhsv(i) 

Else 

PdubO) = dub(i) 

Frvdif 

10 continue 
Call lineq(it) 

Write{6,555) it 

555 Format(lx,’numberof iterations =’i5) 

Write<16,600)it 

C 

Return 

600 format(lx, , number erf solver iterations = \i4) 

End 

Subroutine lineq(it) 
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CJUOUU U U oUU 


Program to solve linear equations (based on: j. Comp. Phisics, 
’The iterative calculation of.../, 17. Pp. 87-94, 1975) 

For information: call charley bauscftlicfrer (415) 694-623 1 

Include paramdat’ 

Include commonJ’ 

Dimension v(ospdim,20), w(nspdim,20),a(20,20),al(20), 

+ buf(nspdim), gg(20), buf 1(400), buf2(nspdim), buf3(nspdim) 


Npan = nspdim 
Solres = 0.000005 
Maxit = 200 
Thresh = solres 
Matdim = npan 
I ms = 0 


Initial guess for starting solution vector 

If(itstep.gL0)then 
Do 10 i=l, matdim 
Buf2(i) = pdub(i) 

10 Continue 
Go to 800 
F.ndif 

Do 20 i=l. matdim 
Ahnn = diag(i) 

If(abs(ahnn).ltl.e-7)ahnn = 1.0e-7 
Buf2(i) s rhsv(i)/ahnn 
20 continue 
800 continue 
Write( 16,600) 

Write( 16,601) 

Write(6,899) 

899 Format( 1 x,'solution iteration historyO 
810 continue 
Ims = ims + ! 

Call normal (buf2, matdim) 

It = it + 1 

Call firmab(buf2,buf3,matdim,iiTu) 

Do 30 i=l, matdim 
W(i,ims) = buf3(i) 

V(i,ims) = bu£2(i) 

30 continue 
Do 40 1, ims 

Do 50 j=l .matdim 
Buf(j) = v(j,i) 

Buf2(j) = w(j,ims) 

50 Continue 

A(i,ims) = sdot(matdim,buf,l,buf2,l) 
lf(i.eq.ims)go to 40 
Do 60 j= 1 .matdim 
Buf(j) = w(j,i) 

Buf2(j) = vij^ms) 

60 Continue 

A(ims,i) = sdot(matdiin,buf, 1 ,buG, 1) 

40 continue 
Do 70 i=l, matdim 
Buf2(i) = v(i,ims) 

70 continue 

Gg(ims) = sdot(matdim 4 tsv, 1 ,buf2, 1) 

Iq = 0 

Do 80 i=l,ims 
Do 90 j=l,ims 
Iq = iq + 1 
Bufl(iq) = a(jj) 
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90 Continue 
Al(i) = gg(i) 

80 continue 

Call gss(bufl ,al,ims,ier) 

IfOer.eq. l)then 
Stop 
Endif 

Cony = abs(al(ims)) 

Ifold = 0 
lf(ims.eq.20)then 
Ifold = ims 

Call zero(buf2,maidim) 

Do 100 i=l,ims 
Do 1 10 j=l,matdim 
Buf(j) = w(j,i) 

1 10 Continue 
Ali = al(i) 

Call saxpy(rmtdin%ali,buf,l,buf2,l) 
100 Continue 

Do 120j=l,matdim 
W(j,l) = buf2(j) 

120 Continue 
X = 0.0 
Y = 0.0 

Do 130 i=l,ims 
Do 140 j=l,ims 
X = x + aG j) * al(i) * al(j) 

140 Continue 

Y = y + gg(i) * al(i) 

1 30 Continue 
A(U) = x 
Gg(l) = y 
Ims = 1 
Endif 

Call zero(buf2,matdim) 

Do 150 i=l,ims 
Do I60j=l,matdim 
Buf(j) = w(j,i) 

160 Continue 
Ali = al(i) 

If(ifoidne.0)ali = 1.0 
Call saxpy(matdim t ali t buf,l,bu£2,l) 
150 continue 
Coax = 0.0 
Ipanel = 1 
Do 170 i=l.matdim 
If(abs(rhsv(i)).lLl.e*7) go to 820 
Q = abs((buf2(i) - rhsv(i))/thsv(i)) 
If(conx.ltq)tben 
Ipanel = I 
Endif 

Coax = amaxl(conx,q) 

820 Continue 

Buf20) = buG(i) - rhsv(i) 

170 continue 

Write( 1 6, 602)it,cony,conx, ipanel 
Noconv = 0 

If(conx.itthresh.and.ifold.eq.O)go to 830 
Noconv ss 1 

If(iLeq.maxit) go to 830 
Do 180 i=l,maldim 
Ahnn = diag(i) 

If(abs(ahnn).lL 1 .0e-7)ahnn = 1 .0e-7 
Buf2(i) = buG(i)/ahnn 
180 continue 

Call normal (buf2,matdim) 

Lp = max0<ifold t ims) 

Do 190 i= Up 



Do 200 j= 1 .matdira 
Buf(j) = v(j,i) 

200 Continue 

X = sdot(matdim,buf, 1 ,buf2, 1) 

Call saxpy(matdim,-x,buf,l,bui2,l) 

Call normal (buf2,matdim) 

190 continue 
If(ifold.eq.O)go to 810 
Call zero(buf3,inatdim) 

Do 210 i=l,ifold 
Do 220 j=l,matdim 
Buf(j) = v(j,i) 

220 Continue 
All = ai(i> 

Call sajtpy(matdim,ali,buf,l,buf3,l) 

210 continue 
Do 230 i=l,matdim 
V(i,l) = bui3(i) 

230 continue 
Al(l)= 1.0 
Go to 810 
830 continue 
Call zero(buf,matdim) 

Do 240 i=l, ims 
Do 250 j=l,matdim 
Buf2(j) = v(j,i) 

250 Continue 
Alii = al(i) 

Call saxpy(matdim,aiil,buf2 i l,buf,l) 

240 continue 
Do 260 i=l ( matdim 
Dub(i) - buf(i) 

260 continue 
lf(noconv.eq.l) then 
WriteO 6,603) 

End if 

600 format(lhl) 

601 format(lx/solution iteration history/) 

602 formate it=\i5/ al(i) 715.8; hv-g 715.8,' panel = \i5) 

603 format(// — no convergence 0 

Return 

End 

Subroutine frmab(a,b,matdim,irawm) 

C 

Include 'paramdat' 

Include 'common^ 

C 

Dimension a(matdim),b<matdim) 

Rewind irawm 
Npan = matdim 
Do 10 i=l,opan 
Ii = I 
B(i) = 0.0 
If(inram.eq.l)then 

Read(ira wm)edubicww(inranvj) j= 1 .npan) 

Ii = l 
Endif 

Do 20 j=l,npan 

BeO = b(i) + a(j) * dubicww(iij) 

20 Continue 
10 continue 
Return 
End 

Subroutine zero(a,len) 

Dimension aOen) 

Do 10 i=l,Ien 
A<i> = 0.0 
10 continue 
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Return 

End 

Subroutine normal(ajen) 

Dimension a(len) 

X = 0.0 
Do 10 i=l,len 
X = x + a(i) * a(i) 

10 continue 
X = 1 .0/sqrt(x) 

Do 20 i=l,len 
A(i) = a(i) * x 
20 continue 
Return 
End 

Subroutine gss(b,g,ninix,ier) 

Dimension b(nmix,nmix),g(nmix) 

Data zerol/1.0e-16/ 

Ier = 0 

Do 10 i=l,nmix 
lf(abs(b(i,i)).lLzerol)go to 800 
Fx = Ub(u) 

Goto 810 
800 Continue 

If(i.eq.nmix)go to 820 
11 = 1+1 
C Pivot section 

Do 20 j=il,nmix 
If(abs(b(j > i)).ltzerol)go to 20 
Fx = Ub(j,i) 

Do 30 l=i,nmix 
Temp = b(j,l) 

B(j,l) = b<U) 

B(i,l) = temp 
30 Continue 
Tmp = g(j) 

G(j) = g(t) 

G(i) = tmp 
Goto810 
20 Continue 
Go to 820 
810 Continue 
G(i) = g(i)*fx 
Do 40 j=i,nmix 
B(i j) = b(ij) * f* 

40 Continue 
Do 50 j=l,nirnx 
lf(i.eq.j)go to 50 
Y = b(j,i) 

G(j) * g(j) - gO) * y 
Do 60 k=i ( nmix 
B(j,k) = tKj.k) - y * b(i,k) 

60 Continue 
50 Continue 
10 continue 
Return 
820 continue 
Write( 16,600) 

600 focmatC Abort invert Singular matrix ) 

Ier= 1 
Return 
End 

Subroutine saxpy(n,sa,sx,incx t sy,incy) 

Constant times a vector plus a vector. 

Uses unrolled loops for increments equal to one. 
Jack dongarra. Unpack, 3/1 1/78. 

Dimension sx(n),sy(n) 
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If(n.le.O)return 
lf(sa.eq.0.0)retum 
If(incx.eq.l.and.incy.eq.l)go to 20 

Code for unequal increments or equal increments 
Not equal to 1 

U=1 

Iy=l 

If(incx.ltO)ix=(-n+ 1 )*incx+ 1 
If(incy.H.0)iy=<-n+ 1 )*incy+ 1 
Do 10 i=l,n 
Sy(iy)=sy(iy)+sa*sx(ix) 

Ix=ix+incx 
Iy=iy+incy 
10 continue 
Return 

Code for both increments equal to 1 


Clean-up loop 

0 m=mod(n,4) 
lf(m.eq.0)go to 40 
Do 30 i=l,m 
Sy(i)=5y(i)+sa*sx(i) 

0 continue 
lf(n.lL4)return 
Ompl=ro+l 
Do 50 i=mp 1,0,4 
Sy(i)=sy(i)+sa*sx(i) 

Sy(i+ l)=sy(i+ l)+sa*sx(i+ 1) 
Sy(i+2)=sy(i+2)+sa*sx(i+2) 
Sy(i+3)=sy(i+3)+sa*sx(i+3) 

0 continue 

Return 

End 

Real function sdot(n^x,incx^y,incy) 

Forms the dot product of two vectors. 

Uses unrolled loops for increments equal to one. 

Jack dongana. Unpack, 3/1 1/78. 

Dimension sx(n)^y(n) 

Stemp=0.0e0 

Sdot=0.0e0 

lf(n.le.0)retuni 

If(incx.eq. 1 .and.incy.eq. I ) go to 20 

Code for unequal increments or equal increments 
Not equal to one. 


Ix= 1 
Iy=l 

If(incxJL0)ix=<-n+l)*incx+l 
[f(incy.lt0)iy=<-Q+ l)*incy+ 1 
Do 10 i=l,n 

Stemp=stemp+ sx(ix) *sy(i y) 
Ix=ix+incx 
Iy=iy+incy 
10 continue 
sdot=stemp 
return 

Code for both increments equal to 1 
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